rm(list = ls(all = TRUE)) #limpiar workspace
setwd("~/Dropbox/REDISTRICTING/KING_MODEL_1991-2012")

# ------------
# Create data
# ------------
vote <- matrix(
  c(0.1682, 0.5847, 0.0791,  #1991
    0.2498, 0.4858, 0.1612,  #1994
    0.2585, 0.3800, 0.2498,  #1997
    0.3824, 0.3692, 0.1868,  #2000
    0.3073, 0.3679, 0.1761,  #2003
    0.3339, 0.2821, 0.2899,  #2006
    0.2801, 0.3832, 0.1220,  #2009
    0.2589, 0.3658, 0.2695), #2012
  ncol=3, byrow=TRUE, 
  dimnames=list(seq(1991,2012,3), c("PAN","PRI","PRD")))

seats <- matrix(
  c(10 , 289,  0,  #1991
    19 , 273,  6,  #1994
    64 , 165, 70,  #1997
    142, 132, 26,  #2000
    80 , 163, 55,  #2003
    137,  65, 98,  #2006
    70 , 188, 39,  #2009
    52 , 174, 71), #2012
  ncol=3, byrow=TRUE, 
  dimnames=list(seq(1991,2012,3), c("PAN","PRI","PRD")))

# Plot seats-votes
pdf("seats_votes.pdf", width=4, height=4)
par(bty="n")
  col <- c(rep("royalblue3", nrow(seats)),
           rep("firebrick", nrow(seats)),
           rep("gold3", nrow(seats)))
  lim <- range(vote, seats/300)
  plot(as.vector(vote), as.vector(seats)/300,
       xlim=lim, ylim=lim,
       main="Elecciones Diputados Federales, \n1991-2012",
       xlab="% Votos", ylab="% Curules",
       pch=16, col=rgb(t(col2rgb(col)), alpha=150, max=255)
  )
  #Labels for election
  text(as.vector(vote), as.vector(seats)/300,
       labels=rep(c("91","94","97","00","03","06","09","12"), ncol(seats)),
       col=col, pos=4, cex=0.3, offset=0.3)
  abline(0,1)
  text(0.7,0.7, expression("" %<-% "Repr. Prop."), pos=4, cex=.7)
dev.off()

#-------------
## Model
#-------------
library(rjags)
library(coda)
library(R2jags)
load.module("glm")
load.module("dic")

set.seed(12345)
fit <- jags(model="01_King_Model.bug",
            data=list(s=seats, v=vote, n=rowSums(seats), 
                      NOBS=nrow(seats), NPARTY=ncol(seats),
                      xstar=0.33),
            n.chains=2, n.iter=15000, n.burnin=5000, n.thin=2,
            parameters.to.save=c("rho", "lambda", "ev", "pv", "Pred"))
#The initial 5,000 MCMC scans were discarded as burn-in. 
#The posterior summaries are based on a posterior sample of size 10,000 
#formed by running each chain for an additional 10,000 scans and storing every 2nd scan.
print(fit)

## Plots in coda
fit.mcmc <- as.mcmc(fit)
pdf("parameters.pdf", width=6, height=8)
plot(fit.mcmc[,c("rho","lambda[2]","lambda[3]")])
dev.off()

## ------------
# Model Fit
## ------------
pdf("fit.pdf", width=4.5, height=4.5)
par(bty="n")
  fitted <- fit$BUGSoutput$median$Pred
  centile <- fit$BUGSoutput$summary
  centile <- centile[grep("Pred", rownames(centile)), c("2.5%","97.5%")]
  #Point Estimates
  plot(as.vector(seats), as.vector(fitted),
       main="Bondad de ajuste",
       xlab="Observed", ylab="Predicted",
       xlim=c(0,300), ylim=c(0,300),
       pch=16, col=rgb(t(col2rgb(col)), alpha=150, max=255)
  )
  #Credible intervals
  segments(as.vector(seats), centile[,1], as.vector(seats), centile[,2],
           col=col, lwd=1.5)
  #Labels for election
  text(as.vector(seats), as.vector(fitted),
       labels=rep(c("91","94","97","00","03","06","09","12"), ncol(seats)),
       col=col, pos=4, cex=0.3, offset=0.3)
  abline(0,1)
dev.off()
